library(leaflet)
popup = c("Robin", "Jakub", "Jannes")
leaflet() %>%
addProviderTiles("Stamen.Toner") %>%
addCircles(lng = c(-3, 23, 11),
lat = c(52, 53, 49),
popup = popup)
#chapter2
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(raster)
## Loading required package: sp
library(spData) # load geographic data
library(spDataLarge) # load larger geographic data
world
## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 11
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melanesia Sove… 1.93e4 8.86e5 70.0
## 2 TZ Tanzania Africa Africa Eastern … Sove… 9.33e5 5.22e7 64.2
## 3 EH Western … Africa Africa Northern… Inde… 9.63e4 NA NA
## 4 CA Canada North Am… Americas Northern… Sove… 1.00e7 3.55e7 82.0
## 5 US United S… North Am… Americas Northern… Coun… 9.51e6 3.19e8 78.8
## 6 KZ Kazakhst… Asia Asia Central … Sove… 2.73e6 1.73e7 71.6
## 7 UZ Uzbekist… Asia Asia Central … Sove… 4.61e5 3.08e7 71.0
## 8 PG Papua Ne… Oceania Oceania Melanesia Sove… 4.65e5 7.76e6 65.2
## 9 ID Indonesia Asia Asia South-Ea… Sove… 1.82e6 2.55e8 68.9
## 10 AR Argentina South Am… Americas South Am… Sove… 2.78e6 4.30e7 76.3
## # … with 167 more rows, and 2 more variables: gdpPercap <dbl>,
## # geom <MULTIPOLYGON [°]>
plot(world)
## Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to plot
## all
names(world)
## [1] "iso_a2" "name_long" "continent" "region_un" "subregion" "type"
## [7] "area_km2" "pop" "lifeExp" "gdpPercap" "geom"
#world["lifeExp"]
world %>% dplyr::select(lifeExp) %>% st_drop_geometry()
## # A tibble: 177 x 1
## lifeExp
## * <dbl>
## 1 70.0
## 2 64.2
## 3 NA
## 4 82.0
## 5 78.8
## 6 71.6
## 7 71.0
## 8 65.2
## 9 68.9
## 10 76.3
## # … with 167 more rows
#st_drop_geometry(world["lifeExp"])
class(world)
## [1] "sf" "tbl_df" "tbl" "data.frame"
plot(world["lifeExp"])
world_asia = world[world$continent == "Asia", ]
asia = st_union(world_asia)
## although coordinates are longitude/latitude, st_union assumes that they are planar
#> although coordinates are longitude/latitude, st_union assumes that they are planar
#world_asia
#asia
#plot(world_asia)
plot(asia)
plot(world["pop"],reset = FALSE)
plot(asia, add = TRUE, col="red")
multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
multi_point <- st_multipoint(multipoint_matrix)
multipoint_matrix
## [,1] [,2]
## [1,] 5 2
## [2,] 1 3
## [3,] 3 4
## [4,] 3 2
class(multi_point)
## [1] "XY" "MULTIPOINT" "sfg"
lnd_point = st_point(c(0.1, 51.5)) # sfg object
class(lnd_point)
## [1] "XY" "POINT" "sfg"
lnd_geom = st_sfc(lnd_point, crs = 4326) # sfc object
class(lnd_geom)
## [1] "sfc_POINT" "sfc"
lnd_attrib = data.frame( # data.frame object
name = "London",
temperature = 25,
date = as.Date("2017-06-21")
)
lnd_attrib
## name temperature date
## 1 London 25 2017-06-21
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom) # sf object
lnd_sf
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 0.1 ymin: 51.5 xmax: 0.1 ymax: 51.5
## Geodetic CRS: WGS 84
## name temperature date geometry
## 1 London 25 2017-06-21 POINT (0.1 51.5)
#load data
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
new_raster = raster(raster_filepath)
new_raster
## class : RasterLayer
## dimensions : 457, 465, 212505 (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333 (x, y)
## extent : -113.2396, -112.8521, 37.13208, 37.51292 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /usr/local/lib/R/site-library/spDataLarge/raster/srtm.tif
## names : srtm
## values : 1024, 2892 (min, max)
plot(new_raster)
multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_raster_file
## [1] "/usr/local/lib/R/site-library/spDataLarge/raster/landsat.tif"
r_brick = brick(multi_raster_file)
r_brick
## class : RasterBrick
## dimensions : 1428, 1128, 1610784, 4 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 301905, 335745, 4111245, 4154085 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs
## source : /usr/local/lib/R/site-library/spDataLarge/raster/landsat.tif
## names : landsat.1, landsat.2, landsat.3, landsat.4
## min values : 7550, 6404, 5678, 5252
## max values : 19071, 22051, 25780, 31961
plot(r_brick)
plotRGB(r_brick[[1:3]]/max(r_brick[[1:3]]),3,2,1,scale=T)
#This is a dataset containing the four bands (2, 3, 4, 5) of the Landsat 8 image for the area of Zion National Park
ndvi = (r_brick[[4]] - r_brick[[3]])/(r_brick[[3]]+r_brick[[4]])
plot(ndvi)
#chapter3
library(sf)
library(raster)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(stringr) # for working with strings (pattern matching)
library(tidyr) # for unite() and separate()
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:raster':
##
## extract
sel_area = world$area_km2 < 10000
small_countries = world[sel_area, ]
small_countries
## Simple feature collection with 7 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -67.24243 ymin: -16.59785 xmax: 167.8449 ymax: 50.12805
## Geodetic CRS: WGS 84
## # A tibble: 7 x 11
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 PR Puerto Ri… North Am… Americas Caribbean Depe… 9225. 3534874 79.4
## 2 PS Palestine Asia Asia Western … Disp… 5037. 4294682 73.1
## 3 VU Vanuatu Oceania Oceania Melanesia Sove… 7490. 258850 71.7
## 4 LU Luxembourg Europe Europe Western … Sove… 2417. 556319 82.2
## 5 <NA> Northern … Asia Asia Western … Sove… 3786. NA NA
## 6 CY Cyprus Asia Asia Western … Sove… 6207. 1152309 80.2
## 7 TT Trinidad … North Am… Americas Caribbean Sove… 7738. 1354493 70.4
## # … with 2 more variables: gdpPercap <dbl>, geom <MULTIPOLYGON [°]>
plot(small_countries["pop"])
world %>% filter(area_km2 < 10000) %>%
dplyr::select(pop) %>%
plot()
##1.GDPperCapitaがtop10の国をピックアップして地図化してみましょう
#top_n
world %>% top_n(n=10,wt=gdpPercap) %>%
dplyr::select(gdpPercap) %>%
plot()
##2.GDPperCapitaがworst10の国をピックアップして地図化してみましょう
#top_n
world %>% top_n(n=-10,wt=gdpPercap) %>%
dplyr::select(gdpPercap) %>%
plot()
##3.日本はGDPperCapitaで上から何番目でしょう
world %>% dplyr::select(name_long,gdpPercap)
## Simple feature collection with 177 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 3
## name_long gdpPercap geom
## <chr> <dbl> <MULTIPOLYGON [°]>
## 1 Fiji 8222. (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135,…
## 2 Tanzania 2402. (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09…
## 3 Western Saha… NA (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27…
## 4 Canada 43079. (((-122.84 49, -122.9742 49.00254, -124.9102 49.9845…
## 5 United States 51922. (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, …
## 6 Kazakhstan 23587. (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48…
## 7 Uzbekistan 5371. (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45…
## 8 Papua New Gu… 3709. (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -…
## 9 Indonesia 10003. (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 …
## 10 Argentina 18798. (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85,…
## # … with 167 more rows
#gdpPercapをソートしたい
#sort(world$gdpPercap,decreasing = T)
gdp_sort <- world %>% arrange(desc(gdpPercap)) %>%
dplyr::select(gdpPercap) %>%
st_drop_geometry()
#日本のgdpPercapを取り出したい
jp_gdp <- world %>% filter(name_long == "Japan") %>%
dplyr::select(gdpPercap) %>%
st_drop_geometry() %>%
as.numeric()
#gdpPercapが上から何番目か調べられれば良い
#which()
which(gdp_sort==jp_gdp)
## [1] 22